This is an open assessment looking at potential health effects of a national fish promotion program in Finland. The details of the assessment are described on Opasnet. This file contains the R code to run the assessment model.
Calculation is based on BONUS GOHERR project and its http://en.opasnet.org/w/Goherr_assessment.
# This is code Op_fi5923/model on page [[Kotimaisen kalan edistämisohjelma]]
# The code was forked from Op_fi5889/model on page [[Ruori]] and Op_en7748/model on page [[Goherr assessment]]
library(OpasnetUtils) # Install the newest version from https://github.com/jtuomist/OpasnetUtils not from CRAN.
library(ggplot2)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
# First empty all objects for a fresh start. Otherwise may be problems with CheckDecisions.
oempty(all=TRUE)
openv.setN(10)
dat <- opbase.data("Op_fi5923", subset="Malliparametrit")[-1]
dec <- opbase.data("Op_fi5923", subset="Decisions")[-1]
DecisionTableParser(dec)
CTable <- opbase.data("Op_fi5923",subset="CollapseMarginals")
for(i in 1:ncol(CTable)) {CTable[[i]] <- as.character(CTable[[i]])}
CollapseTableParser(CTable)
objects.latest("Op_en2031", code_name="subgrouping") # [[Exposure-response function]] subgrouping
## Loading objects:
## subgrouping
cat("Laskennassa käytetty data.\n")
## Laskennassa käytetty data.
dat
cat("Tarkastellut päätökset.\n")
## Tarkastellut päätökset.
dec
cat("Aggregoidut marginaalit.\n")
## Aggregoidut marginaalit.
CTable
#' prepare adjusts the data table for ovariables. Requires function subgrouping from code Op_en2031/initiate on page [[Exposure-response function]]
#' @param dat data.frame
#' @param type type of data that is used. Must match content in column Type
#' @param drop columns to remove
#' @return data.frame
prepare <- function(dat, type=NULL, drop=NULL) {
out <- dat
if(!is.null(type)) out <- out[out$Type %in% type , ]
if(!is.null(drop)) out <- out[!colnames(out) %in% drop]
return(subgrouping(out))
}
amount <- Ovariable("amount", data = prepare(dat, "amount", c("Type","Response","Exposure_agent"))) # Filleted weight, i.e. no loss.
# Exposure:To child and To eater not needed, because dioxins are not (yet) included
conc_vit <- Ovariable(
"conc_vit",
ddata = "Op_en1838",
subset = "Fineli data for common fish species"
)
# dependencies = data.frame(Name="dummy"),
# formula = function(...) {
# df = opbase.data(
# "Op_en1838", # [[Concentrations of beneficial nutrients in fish#Fineli]]
# subset = "Fineli data for common fish species"
# )
df = conc_vit@data
df$Nutrient[df$Nutrient=="D-vitamiini (µg)"] <- "Vitamin D"
df$Nutrient[df$Nutrient=="rasvahapot n-3 moni-tyydyttymättömät (g)"] <- "Omega3"
df$Nutrient[df$Nutrient=="rasvahappo 18:3 n-3 (alfalinoleenihappo) (mg)"] <- "ALA"
df$Nutrient[df$Nutrient=="rasvahappo 22:6 n-3 (DHA) (mg)"] <- "DHA"
df$Nutrient[df$Nutrient=="proteiini (g)"] <- "Fish"
df$conc_vitResult[df$Nutrient=="Fish"] <- "1"
df <- dropall(df[df$Nutrient %in% c("Vitamin D", "Omega3", "ALA", "DHA", "Fish") , ])
# return(df)
# }
#)
conc_vit@data <- df
exposure <- Ovariable(
"exposure",
dependencies = data.frame(Name=c("conc_vit", "amount")),
formula = function(...) {
# First, match KKE-classification with Fineli classification
conc_vit$Kala <- NULL
tmp <- Ovariable(
output = data.frame(
Kala = c("Kasvatettu", "Kaupallinen", "Kirjolohi", "Silakka", "Vapaa-ajan", "Muu tuonti", "Tuontikirjolohi", "Tuontilohi"),
Fish = c("Whitefish", "Average fish","Rainbow trout", "Herring", "Average fish", "Average fish", "Rainbow trout", "Salmon"),
Result = 1
))
amount <- amount * 1000 / 5.52 /365.25 # M kg/a per 5.52M population --> g/d per person.
out <- conc_vit * tmp * amount
colnames(out@output)[colnames(out@output)=="Nutrient"] <- "Exposure_agent"
return(out)
},
marginal = c(TRUE, TRUE, FALSE)
)
# Incidence-based data not needed
frexposed <- 1
population <- 0
incidence <- 0
ERFChoice <- Ovariable(
"ERFchoice",
data=data.frame(
Response=c("Loss in child's IQ points","CHD2 mortality","Breast cancer","All-cause mortality","Depression"),
Exposure_agent=c("DHA","Omega3","Omega3","Fish","Fish"),
Result=1)
)
case_burden <- Ovariable("case_burden", data= prepare(dat,"case burden",c("Type","Exposure_agent","Unit")))
objects.latest("Op_en5917", code_name="InpBoD") # [[Disease risk]] InpBoD for GOHERR
## Loading objects:
## InpBoD
InpBoD <- EvalOutput(InpBoD) # Evaluated because is not a dependency but an Input
InpBoD <- InpBoD[InpBoD$Country=="FI",colnames(InpBoD@output)!="Country"]
objects.latest("Op_en2261",code_name="BoDattr2") # [[Health impact assessment]]
## Loading objects:
## BoDattr
tryCatch(BoDattr <- EvalOutput(BoDattr, verbose=TRUE))
## Evaluating BoDattr ...
##
## - BoD fetched successfully!
##
## - PAF fetched successfully!
## - Evaluating BoD ...
## - - Evaluating case_burden ...
##
## done(0.02 secs)!
## - - Checking case_burden marginals ... Response, Iter, case_burdenSource recognized as marginal(s).
## Loading required package: reshape2
## - - Processing case_burden marginal collapses ... done!
##
## - done(0.12 secs)!
## - Checking BoD marginals ... Response, Iter, BoDSource recognized as marginal(s).
## - Processing BoD inputs ... done!
## - Processing BoD marginal collapses ...
## Warning in oapply(variable, FUN = fun[[i]], cols = cols[[i]], na.rm = TRUE):
## While oapplying BoD, found NAs in indices: Group, InpBoDSource. They were
## automatically filled using fillna, which may result in a multiplied population.
## Please check your ovariable before using oapply.
## done!
## - Evaluating PAF ...
##
## - - dose fetched successfully!
##
## - - ERF fetched successfully!
##
## - - RR fetched successfully!
##
## - - P_illness fetched successfully!
##
## - - sumExposcen fetched successfully!
##
## - - mc2d fetched successfully!
##
## - - mc2dparam fetched successfully!
## - - Evaluating dose ...
##
## - - - BW fetched successfully!
## - - - Evaluating exposure ...
## - - - - Evaluating conc_vit ...
##
## done(0 secs)!
## - - - - Checking conc_vit marginals ... Kala, Fish, Nutrient, Iter, conc_vitSource recognized as marginal(s).
## - - - - Processing conc_vit decisions ... done!
## - - - - Evaluating amount ...
##
## done(0 secs)!
## - - - - Checking amount marginals ... Kala, Scenario, amountSource recognized as marginal(s).
##
## --- done(0.26 secs)!
## - - - Checking exposure marginals ... Fish, Exposure_agent, Iter, conc_vitSource, Adjust, Scenario, amountSource, exposureSource recognized as marginal(s).
## - - - Processing exposure marginal collapses ... done!
## - - - Evaluating BW ...
##
## done(0 secs)!
## - - - Checking BW marginals ... BWSource recognized as marginal(s).
##
## -- done(25.94 secs)!
## - - Checking dose marginals ... Exposure_agent, Iter, conc_vitSource, Adjust, Scenario, amountSource, Scaling, exposureSource, BWSource, doseSource recognized as marginal(s).
## - - Processing dose marginal collapses ... done!
## - - Evaluating ERF ...
##
## - - - ERF_env fetched successfully!
##
## - - - ERF_omega3 fetched successfully!
##
## - - - ERF_mehg fetched successfully!
##
## - - - ERF_diox fetched successfully!
##
## - - - ERF_vit fetched successfully!
##
## - - - ERF_micr fetched successfully!
##
## - - - ERFchoice fetched successfully!
## - - - Evaluating ERF_env ...
##
## done(0.01 secs)!
## - - - Checking ERF_env marginals ... Exposure_agent, Response, Subgroup, Exposure, ER_function, Scaling, Exposure_unit, Observation, Iter, ERF_envSource recognized as marginal(s).
## - - - Evaluating ERF_omega3 ...
##
## done(0.01 secs)!
## - - - Checking ERF_omega3 marginals ... Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling, Observation, Iter, ERF_omega3Source recognized as marginal(s).
## - - - Evaluating ERF_mehg ...
##
## done(0 secs)!
## - - - Checking ERF_mehg marginals ... Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling, Observation, Iter, ERF_mehgSource recognized as marginal(s).
## - - - Evaluating ERF_diox ...
##
## done(0.01 secs)!
## - - - Checking ERF_diox marginals ... Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling, Observation, Iter, ERF_dioxSource recognized as marginal(s).
## - - - Evaluating ERF_vit ...
##
## done(0 secs)!
## - - - Checking ERF_vit marginals ... Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling, Observation, ERF_vitSource recognized as marginal(s).
## - - - Evaluating ERF_micr ...
##
## done(0 secs)!
## - - - Checking ERF_micr marginals ... Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling, Observation, ERF_micrSource recognized as marginal(s).
## - - - Evaluating ERFchoice ...
## - - - - Processing exposure marginal collapses ... done!
##
## --- done(0.05 secs)!
## - - - Checking ERFchoice marginals ... Exposure_agent, ERFchoiceSource recognized as marginal(s).
##
## -- done(2.86 mins)!
## - - Checking ERF marginals ... Exposure_agent, Response, Exposure, ER_function, Scaling, Observation, Iter, ERFSource recognized as marginal(s).
## - - Processing ERF marginal collapses ... done!
## - - Evaluating RR ...
## - - - Processing dose marginal collapses ... done!
## - - - Processing ERF marginal collapses ... done!
##
## -- done(0.24 secs)!
## - - Checking RR marginals ... Exposure_agent, Response, ER_function, Scaling, Iter, ERFSource, conc_vitSource, Adjust, Scenario, amountSource, doseSource, RRSource recognized as marginal(s).
## - - Evaluating P_illness ...
##
## done(0 secs)!
## - - Checking P_illness marginals ... Response, Illness, Age, P_illnessSource recognized as marginal(s).
##
## - done(6.1 mins)!
## - Checking PAF marginals ... Exposure_agent, Response, ER_function, Scaling, Iter, ERFSource, conc_vitSource, Adjust, Scenario, amountSource, doseSource, RRSource, PAFSource recognized as marginal(s).
## - Processing PAF marginal collapses ... done!
##
## done(6.93 mins)!
## Checking BoDattr marginals ... Response, Iter, Group, InpBoDSource, Exposure_agent, conc_vitSource, Adjust, Scenario, amountSource, PAFSource, BoDattrSource recognized as marginal(s).
###################
# Graphs
trim <- function(ova) return(oapply(ova, NULL, mean, "Iter")@output)
plot_ly(trim(amount), x=~Scenario, y=~amountResult, color=~Kala, type="bar") %>%
layout(yaxis=list(title="Kalan kokonaiskulutus Suomessa (milj kg /a)"), barmode="stack")
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
plot_ly(trim(conc_vit), x=~Nutrient, y=~conc_vitResult, color=~Kala, type="scatter", mode="markers") %>%
layout(yaxis=list(title="Concentrations of nutrients (mg or ug /g)"))
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
tmp <- exposure / Ovariable(
output = data.frame(
Exposure_agent = c("Fish","Vitamin D", "Omega3", "ALA", "DHA"),
Result = c(1, 1, 1000, 1000, 1000)
), marginal = c(TRUE, FALSE)
)
plot_ly(trim(tmp), x=~Scenario, y=~Result, color=~Exposure_agent, text=~Exposure_agent, type="bar") %>%
layout(yaxis=list(title="Exposure to nutrients (g or ug /d)"))
summary(conc_vit, "mean")
summary(exposure,"mean")
summary(amount,"mean")
summary(BoD,"mean")
summary(BoDattr,"mean")
cat("Kalaperäisiä tautitaakkoja Suomessa\n")
## Kalaperäisiä tautitaakkoja Suomessa
tmp <- summary(oapply(BoDattr,NULL,sum,"Group")[BoDattr$Scenario=="BAU",])
data.frame(
Altiste = tmp$Exposure_agent,
Vaikutus = tmp$Response,
Keskiarvo = as.character(signif(tmp$mean,2)),
"95 luottamusväli" = paste0(signif(tmp$Q0.025,2)," - ", signif(tmp$Q0.975,2)),
Keskihajonta = signif(tmp$sd,2)
)#[rev(match(lev, tmp$Exposure_agent)),]
ggplot(trim(BoDattr), aes(x=Scenario, weight=BoDattrResult, fill=Response))+geom_bar()
ggplot(trim(BoDattr), aes(x=Scenario, weight=BoDattrResult, fill=Exposure_agent))+geom_bar()
plot_ly(trim(BoDattr), x=~Scenario, y=~BoDattrResult, color=~Response, text=~paste(Group, Exposure_agent, sep=": "), type="bar") %>%
layout(yaxis=list(title="Disease burden (DALY /a); CHD2=coronary heart disease"), barmode="stack")
## Warning: Ignoring 12 observations
## Warning: n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning: n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
tmp <- summary(oapply(BoDattr,NULL,sum,c("Age","Exposure_agent")))
data.frame(
Altiste = tmp$Response,
Keskiarvo = signif(tmp$mean,2),
"95 luottamusväli" = paste0(signif(tmp$Q0.025,2)," - ", signif(tmp$Q0.975,2)),
Keskihajonta = signif(tmp$sd,2)
)
################ Insight network
gr <- scrape(type="assessment")
objects.latest("Op_en3861", "makeGraph") # [[Insight network]]
## Loading objects:
## makeGraph
gr <- makeGraph(gr)
## Loading required package: DiagrammeR
## Loading objects:
## formatted
## Loading objects:
## chooseGr
#export_graph(gr, "ruori.svg")
render_graph(gr)
##################### Diagnostics
objects.latest("Op_en6007", code_name="diagnostics")
## Loading objects:
## showind
## binoptest
## showLoctable
## ovashapetest
showLoctable()
showind()
## incidence is not an ovariable.
## population is not an ovariable.
## incidence is not an ovariable.
## frexposed is not an ovariable.
## sumExposcen is not an ovariable.
## mc2d is not an ovariable.
## mc2dparam is not an ovariable.